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Abstract 

This  work  considers  the  problem  of  estimating  directions  of  arrival  (DOAs)  from  sensor-array 
data  when  the  positions  of  the  sensors  are  not  precisely  known.  The  proposed  algorithm  is 
an  extension  of  a  previously  developed  approach  based  on  weighted  subspace  fitting.  The  new 
algorithm  shows  improved  robustness  with  respect  to  sensor  position  perturbations,  partic¬ 
ularly  when  estimating  the  DOAs  of  weak  sources.  The  robust  algorithm  treats  the  sensor 
position  perturbations  as  unknown  deterministic  parameters  and  estimates  them  jointly  with 
the  DOAs.  The  distances  between  each  sensor  in  the  perturbed  array  were  constrained  to 
equal  their  nominal  values.  The  algorithm  adds  one  source  at  a  time  to  the  model  until  the 
amplitude  of  an  identified  source  is  below  a  threshold.  For  each  source  added  to  the  model,  the 
algorithm  performs  a  one-dimensional  grid  searches  followed  by  multidimensional  optimization 
over  a  small  region  of  parameter  space.  This  method  is  effective  for  estimating  a  large  number 
of  parameters. 


1  Introduction 


The  performance  of  array  processing  algorithms  is  degraded  when  the  sensor  positions  are 
perturbed  from  their  nominal  values.  Various  methods  have  been  proposed  to  deal  with  this 
problem.  Gilbert  and  Morgan  [1]  considered  the  calculation  of  beamforming  weights.  They 
showed  that  weights  required  to  produce  a  very  narrow  main  beam  gave  a  severely  degraded 
beampattern  when  the  sensor  positions  were  perturbed.  They  derived  a  formula  to  calculate 
the  narrowest  possible  main  beam  subject  to  a  constraint  on  the  degradation  caused  by  sensor 
position  perturbation. 

Schultheiss  and  Ianiello  [2]  considered  the  effect  of  sensor  postion  perturbations  on  the 
problem  of  estimating  the  range  and  bearing  of  a  single  source.  They  showed  that,  to  first 
order,  conventional  beamforming,  which  is  optimal  when  there  is  no  perturbation,  remains 
optimal  for  small  array  perturbations.  However,  this  result  is  valid  only  for  a  single  source. 

In  this  paper  we  propose  a  robust  algorithm  for  estimating  the  directions  of  arrival  of 
multiple  sources  based  on  weighted  subspace  fitting.  This  is  an  extension  of  the  WSF  algorithm 
presented  in  [3].  In  addition  to  presenting  a  cost  function,  [3]  showed  how  to  obtain  initial 
DOA  estimates,  determine  the  number  of  sources  in  the  data,  and  handle  multiple  clusters  of 
closely  spaced,  non-resol vable  sources.  However,  the  algorithm  in  [3]  assumed  that  the  sensor 
positions  were  precisely  known,  and  its  performance  is  degraded  when  there  is  sensor  position 
perturbation.  In  this  paper,  we  remove  to  some  degree  this  degradation  by  jointly  estimating 
the  post  ion  perturbations,  which  are  considered  to  be  deterministic  unknown  parameters,  and 
the  DOAs. 

There  is  related  work  in  references  [4]- [6].  In  [4|  a  WSF  cost  function  was  derived  based 
on  the  asymptotic  statistics  of  perturbed  eigenvectors,  whereas  the  WSF  cost  function  used  in 
this  paper  was  derived  from  a  subspace  perturbation  expansion  [3] .  References  [5]  and  [6]  use 
statistical  models  for  array  perturbations,  including  sensor  position  errors  and  sensor  gain  and 
phase  errors.  A  WSF  cost  function  is  derived  in  which  the  weights  are  calculated  to  minimize 
the  DOA  estimation  error  variance  with  respect  to  both  the  additive  noise  in  the  data  and  the 
array  perturbations. 

2  Data  Model 

The  model  for  the  noise-free  signal  at  a  single  frequency  is 

Y  =  A(0o,x,y,z)S  =  [ax  ar]S  (1) 

where  Y  ismxn,  A  is  m  x  and  S  is  rx  n.  In  this  application,  m  is  the  number  of  sensors,  n 
is  the  number  of  snapshots  of  array  data,  and  r  is  the  number  of  narrowband  signals  impinging 
on  the  array.  The  vector  of  possible  DOAs  is  6  and  6q  denotes  the  actual  DOAs.  The  vectors 
x,  y,  and  z  contain  the  coordinates  (positions)  of  the  sensors. 

In  this  paper  we  consider  a  nominally  uniform  linear  array  of  sensors  whose  coordinates  are: 
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Thus,  the  sensors  lie  along  the  2  axis  with  a  spacing  of  d  =  1/2  between  each  sensor.  The 
units  of  d  are  wavelengths  corresponding  to  an  array  design  frequency  /q. 

The  columns  of  A  are  the  array  manifold  vectors,  and  they  take  the  following  functional 
form: 

a.  _  g—j2ira(z  cos  sin  0*)  ^ 

where  a  is  the  signal  frequency  expressed  as  a  fraction  of  /o,  the  array  design  frequency.  In 
other  words,  the  actual  signal  frequency  is  afo.  Note  that  for  simplicity,  all  of  the  analysis  in 
this  paper  takes  place  in  the  xz  plane  (all  y  components  equal  zero).  The  approach  can  be 
extended  to  a  full  three-dimensional  analysis. 

We  rewrite  the  array  manifold  vectors  in  terms  of  two  vectors 
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that  represent  perturbations  in  the  2  and  x  components  of  sensor  position,  respectively.  We 
assume  that  all  sensor  positions  are  referenced  to  the  first  sensor  (i.e.  there  is  no  perturbation 
associated  with  the  first  sensor).  Thus,  the  actual  sensor  locations  are  as  follows: 
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The  nominal  and  perturbed  sensor  locations  for  the  first  three  sensors  of  a  line  array  are  shown 
in  Fig.  1.  We  treat  the  position  perturbations  in  the  x  direction  (perpendicular  to  the  line  of 
the  array)  as  free  variables.  The  perturbations  in  the  z  direction  (along  the  line  of  the  array) 
are  calculated  so  as  to  keep  the  spacing  between  the  sensors  equal  to  d.  This  is  expressed  with 
the  following  constraint  equations: 

(5xi  —  8xi~  1)2  +  (d  +  Szi  —  dzi-1)2  =  d2,  i  =  1,  ■  ■  ■ ,  rrt  —  1.  (6) 


These  equations  can  be  solved  for  5z  in  terms  of  5x,  with  5x o  =  Szq  —  0,  as  follows 

5zi  =  5z{- 1  +  y/d?  -  ( Sxi  -  Sxi-i)2  -  d,  i  =  (7) 


The  position  vectors  for  the  actual  sensor  locations  are 


x  =  5xa,  y  =  0,  z  =  z0  +  Sza.  (8) 

Thus,  the  array  manifold  vectors  can  be  written  in  terms  of  the  position  perturbation  vectors 
as  follows: 


The  SVD  of  Y  is 
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x  =  nominal  sensor  location 

o  =  actual  sensor  location 

dx  is  a  uniformly  distributed  random  variable 

dz  is  calculated  to  keep  the  sensor-to-sensor  distance  fixed 

Figure  1:  Nominal  and  perturbed  sensor  locations  for  a  3-element  array. 


where  Ui  has  r  columns.  The  columns  of  Ui  and  A(6q)  span  the  same  subspace,  and  therefore 
columns  of  U2  are  orthogonal  to  columns  of  A(6q). 

The  observed  (noisy)  data  is 

Y  =  Y  +  N  (11) 

where  the  elements  of  N  are  taken  to  be  zero  mean  i.i.d.  complex  Gaussian  random  variables 
with  variance  a2  (real  and  imaginary  parts  are  uncorrelated).  The  SVD  of  Y  is 
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where  Ui  has  r  columns. 

The  subspace-fitting  criterion  used  here  is  based  on  the  fact  that,  in  the  noise-free  case, 

U2U^  A(@0,  6xa,  0,  z0  +  Sza )  =  0. 


With  noisy  data  the  previous  expression  will  not  equal  zero  and  we  look  for  the  parameter 
vector  0  that  minimizes  the  equation  error: 


6  =  argmin 
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. 


(12) 


where  (7)  is  used  to  calculate  Sz  as  a  function  of  5x,  and  the  norm  is  defined  as 

||  •  ||^  =  vec(-)"W  vec(-) 
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and  W  is  a  weight  matrix. 

The  weight  matrix  W  is  calculated  from  the  statistics  of  the  vectorized  cost  function  eval¬ 
uated  at  the  true  parameters  [3].  Note  that  the  perturbation  of  the  cost  function  at  the  true 
parameters  is  due  only  to  the  additive  noise.  Thus,  the  matrix  W  in  this  paper  is  the  same  as 
that  used  in  [3].  After  substituting  in  the  weight  matrix,  the  cost  function  may  be  simplified 
to 

C(0,Sx)  =  ||UfA(0)(Uf  A(6»))-1S1||^,  (13) 

where  A  —  A(0,<5x,O,zo  +  <5z),  Ei  =  (Sj  -  <72I)°‘5,  and  <j2  is  the  average  of  the  squared 
singular  values  in  £2-  Finally,  the  subscript  (F5  in  (13)  denotes  the  Frobenious  matrix  norm. 

3  Review  of  Previous  Algorithm 

The  steps  of  the  original,  nonrobust  algorithm  [3]  are  reviewed  in  this  section.  The  modifi¬ 
cations  which  are  needed  to  give  robustness  to  sensor  position  perturbations  are  given  in  the 
following  section.  In  the  original  algorithm,  the  cost  function  is  given  by  (13)  with  the  position 
perturbation  parameters  fixed  at  zero.  The  algorithm  proceeds  as  follows: 

Step  1 

Set  r  —  1  (look  for  one  source)  and  plot  the  reciprocal  of  the  cost  function  on  a  grid  of  points 
to  find  the  maximum.  The  angle  6  that  maximizes  the  reciprocal  of  the  cost  function  is  called 
01. 

Step  2 

To  search  for  the  second  source,  set  r  =  2  and 


where  6\  is  fixed  at  its  value  from  the  first  step.  The  reciprocal  of  the  cost  function  is  plotted 
as  a  function  of  6*2-  The  value  of  60  that  maximizes  the  reciprocal  of  the  cost  function  is  called 
02- 

Step  3 
With 

•-[£. 

as  an  initialization,  C(9)  is  minimized  (e.g.  using  the  Matlab  Optimization  Toolbox),  and  the 
signal  powers  and  noise  powers  are  estimated  [3]. 

Step  4 

This  process  of  adding  one  source  at  a  time  is  repeated  until  the  SNR  estimated  for  any  source 
is  less  than  a  user-specified  threshold  SNR.  The  initial  value  for  each  new  source  is  obtained 
using  a  one-dimensional  grid  search  of  the  cost  function  with  all  previously  found  DO  As  fixed. 
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After  the  initial  value  of  the  new  source  is  found,  the  cost  function  is  minimized  with  respect  to 
all  of  the  DOAs.  This  minimization  is  done  with  tight  upper  and  lower  bounds  on  the  DOAs, 
centered  on  their  current  values. 

Because  of  the  weighting  in  the  cost  function,  the  sources  are  found  roughly  in  the  order  of 
their  powders,  from  strongest  to  weakest.  This  ordering  is  occasionally  altered  when  there  are 
sources  of  nearly  equal  power. 

4  Robust  Algorithm 

Simulation  results  show  that  the  original  algorithm  has  some  robustness  to  sensor  perturba¬ 
tions.  The  main  effects  of  these  perturbations  are  that  the  variances  of  the  DOA  estimates 
increase  and  the  estimated  SNRs  decrease  due  to  the  mismatch  between  the  array  manifold 
model  used  to  calculate  the  SNR  and  the  actual  (perturbed)  array  manifold.  The  original 
algorithm  gives  reasonable  estimates  for  all  but  the  weakest  sources.  Thus,  we  can  use  the 
original  algorithm  with  a  high  SNR  threshold  to  initialize  the  robust  algorithm. 

For  example,  in  the  simulations  shown  in  the  next  section,  the  original  algorithm  uses  an 
SNR  threshold  of  -10  dB  when  there  are  no  sensor  position  perturbations.  In  order  to  deal 
with  sensor  perturbations,  the  threshold  is  set  to  10  dB,  After  sources  above  10  dB  have  been 
estimated,  the  remaining  sources  are  estimated  jointly  with  the  perturbation  parameters  Sx. 
The  grid  searches  for  each  new  DOA  are  done  with  the  Sx  parameters  fixed  at  their  current 
values.  The  minimizations  are  done  with  respect  to  all  DOA  and  perturbation  parameters. 

The  robust  algorithm  treats  the  actual  sensor  position  perturbations  as  unknown  determin¬ 
istic  parameters.  The  initial  guesses  for  the  perturbation  parameters  are  all  zero,  and  tight 
upper  and  lower  bounds  (±0.1)  are  used  for  each  element  of  Sx.  The  perturbation  parameters 
are  expressed  as  a  fraction  of  the  sensor  spacing,  d .  We  use  the  same  sequential  precedure  as 
before,  adding  one  source  at  a  time. 

5  Simulation  Example 

In  this  section  we  consider  a  challenging  simulation  example  consisting  of  seven  moving  sources, 
some  strong  and  others  weak,  at  normalized  frequency  0.4.  The  array  is  a  48-element  uniform 
lineai’  array. 

Fig.  2  shows  a  time-bearing  plot  of  the  estimated  sources  obtained  with  the  original  algo¬ 
rithm;  the  data  are  generated  with  no  array  perturbation.  The  colorbar  gives  the  estimated 
signal-to-noise  ratios  of  each  estimated  source.  The  given  data  contained  1800  snapshots  from 
the  array  and  these  were  processed  in  blocks  of  15  snapshots.  In  this  example  the  number  of 
sources  was  not  assumed  to  be  known  in  advance.  For  each  matrix  processed,  sources  were 
estimated  one  at  a  time.  Each  time  a  new  source  was  estimated,  the  SNR  associated  with  all 
of  the  sources  was  also  estimated.  Additional  sources  were  added  to  the  model  until  the  SNR 
associated  with  one  of  the  estimated  sources  was  less  than  -10  dB. 

In  the  remaining  examples,  the  data  were  generated  with  some  amount  of  array  perturba¬ 
tion.  The  Sx  parameters  specifying  the  perturbed  array  were  generated  as  iid  random  variables 
uniformly  distributed  on  the  interval  [—A,  A],  where  A  is  specified  as  a  fraction  of  the  inter¬ 
sensor  spacing,  d.  Fig.  3  shows  the  result  of  the  original  algorithm  processing  data,  generated 
with  A  =  0.02.  Notice  that  the  weak  sources  near  0  dB  are  not  well  estimated  (compare  with 
Fig.  2). 
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Figure  2:  Time-bearing  plots  for  a  seven-source  simulation.  Color  bar  indicates  estimated 
signal-to- noise  ratio  in  dB.  Original  algorithm.  No  array  perturbation. 


Original  Algorithm,  Deita=0.02 


Figure  3:  Time-bearing  plot  for  seven-source  example.  Original  algorithm.  Data  generated 
with  array  perturbation,  A  —  0.02. 
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Fig.  4  shows  the  result  of  the  original  algorithm  processing  data  generated  with  A  -  0.04. 
The  weak  sources  are  no  longer  visible. 


Original  Algorithm,  Delta=0.04 
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Figure  4:  Time-bearing  plot  for  seven-source  example.  Original  algorithm.  Data  generated 
with  array  perturbation,  A  =  0.04. 

Figs.  5  and  6  show  the  results  of  the  robust  algorithm  processing  data  generated  with 
A  =  0.02  and  0.04,  respectively.  Although  there  is  some  degradation  in  the  estimation  of  the 
weak  sources,  there  is  a  noticeable  improvement  over  the  original  algorithm  (Figs.  3  and  4). 

6  Summary 

The  sensitivity  of  the  WSF  algorithm  described  in  [3]  to  sensor  position  perturbations  can  be 
reduced  by  estimating  these  perturbations  along  with  the  DOAs.  Using  this  approach,  one  is 
able  to  incorporate  constraints  on  the  perturbations.  In  this  paper  the  distances  between  each 
sensor  in  the  perturbed  array  were  constrained  to  equal  their  nominal  values.  The  sequential, 
add  one  source  at  a  time,  procedure  of  one-dimensional  grid  searches  followed  by  multidimen¬ 
sional  optimization  over  a  small  region  of  parameter  space  is  an  effective  method  for  estimating 
a  large  number  of  parameters. 
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Robust  Algorithm,  Delta=0.Q2 
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Figure  5:  Time-bearing  plot  for  seven-source  example.  Robust  algorithm.  Data  generated 
with  array  perturbation,  A  =  0.02. 


Robust  Algorithm,  Delta=0.04 


Figure  6:  Time-bearing  plot  for  seven-source  example.  Robust  algorithm.  Data  generated 
with  array  perturbation,  A  =  0.04. 
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